function fig1(filenumber)
count=filenumber;

N=readdata(count,1);
S=readdata(count,2);
T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
time=readdata(count,6);
%     dtnom=readdata(count,7);
dx=readdata(count,8);
dz=readdata(count,9);
mx=readdata(count,10);
mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
mtype=readdata(count,14);
%mTemp=readdata(count,15);
%mstressxx=readdata(count,16);
%mstressxz=readdata(count,17);
%mstrainxx=readdata(count,18);
%mstrainxz=readdata(count,19);
%mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
%P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

time/(60*60*365*24*1e6)


x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end
figure(1); clf;

mTemp=readdata(count,15);
[grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,2,2);
avtemp=0*zgrid;
for t=1:length(avtemp);
    avtemp(t)=mean(grid(:,t));
end
load colormapfig1a;
ax(1)=subplot(5,3,1:3);%axes('pos',[0.04 0.54 0.42 0.42]);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); ind1=find(Zgrid<40e3); grid(ind1)=nan; 
pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(gca,cmapfig1a); colorbar; axis([0 2000 -40 670]); title(['Temperature (^{\circ}C)']); set(gca,'xticklabel',[]);  ylabel('Depth (km)');
freezeColors; freezeColors(colorbar);

for s=1:length(zgrid)
    grid(s,:)=grid(s,:)-avtemp';
end
ax(2)=subplot(5,3,4:6) ;pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(gca,'jet'); colorbar; axis([0 2000 -40 670]); title(['Deviation from horisontally averaged temperature (^{\circ}C)']); set(gca,'xticklabel',[]);  ylabel('Depth (km)');
freezeColors; freezeColors(colorbar);

mk=zeros(1,N);
mk(find(mtype==1))=2.5; mk(find(mtype==2))=2.5; mk(find(mtype==3 | mtype==4))=4.08*power((298./(mTemp(find(mtype==3 | mtype==4))+273)),0.406); mk(find(mtype==5))=2.5;
[gridk,xgrid,zgrid]=bilingrid(mk,mx,mz,dx,dz,2,2);
[gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,2,2);
gridq=zeros(length(xgrid),length(zgrid)); dxgrid=diff(xgrid); dzgrid=diff(zgrid);
for s=1:length(xgrid)
    for t=2:length(zgrid)-1
        gridq(s,t)=gridk(s,t)*(gridtemp(s,t+1)-gridtemp(s,t-1))/(dzgrid(t-1)+dzgrid(t));
    end
end
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); ind1=find(Zgrid<40e3); gridq(ind1)=nan; 
load colormapfig1c;
ax(3)=subplot(5,3,7:9); pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,gridq); shading flat; set(gca,'ydir','rev'); axis equal; caxis([-0.02 0.07]); colormap(gca,cmapfig1c); colorbar; axis([0 2000 -40 670]); title('Vertical conductive heat flux (W/m^2)'); set(gca,'xticklabel',[]);  ylabel('Depth (km)');
freezeColors; freezeColors(colorbar);

mstressxx=readdata(count,16);
mstressxz=readdata(count,17);
mstrainxx=readdata(count,18);
mstrainxz=readdata(count,19);

load colormapfig1d;
[grid,xgrid,zgrid]=bilingrid(log10(0.5*(mstressxx.^2+mstressxz.^2).^0.5./(mstrainxx.^2+mstrainxz.^2).^0.5),mx,mz,dx,dz,2,2); clear('mstrainxx','mstrainxz','mstressxx','mstressxz');
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; 
ax(4)=subplot(5,3,10:12); pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(gca,cmapfig1d); colorbar; axis([0 2000 -40 670]); title('log Effective viscosity (log(Pa s))'); xlabel('Distance (km)'); ylabel('Depth (km)');
freezeColors; freezeColors(colorbar);

% subplot(5,1,5); pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(gca,'jet'); colorbar; axis([0 2000 -40 670]); title('Strain rate'); xlabel('Distance (km)'); ylabel('Depth (km)');
% pos=get(gca,'pos'); delete(gca); axes('pos',[0.04 0.05 0.42 0.42]);
ax(5)=subplot(5,3,14);
plot(adiabat(zgrid,40e3),(zgrid-40e3)/1e3,'r','linewidth',2); hold on; set(gca,'ydir','rev');
plot(avtemp,(zgrid-40e3)/1e3,'linewidth',2); axis([0 1900 0 670]); xlabel('Temperature(^{\circ}C)'); ylabel('Depth (km)');
a=legend('Theoretical isentropic temperature at \newline potential temperature T_p=1315 ^{\circ}C','Horisontally averaged temperature ','location','southwest');
set(a,'fontsize',6);

